# Jointly testing summer vs. winter (for CARE and non-CARE)
R = matrix(
  data = c(
    -1, 1, 0, 0, rep(0, 4),
    0, 0, -1, 1, rep(0, 4)
  ), nrow = 2, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 2, ncol = 1)
K = both_city[[5]]$nparams
n = both_city[[5]]$nobs
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Jointly testing CARE vs. non-CARE (for each season)
R = matrix(
  data = c(
    -1, 0, 1, 0, rep(0, 4),
    0, -1, 0, 1, rep(0, 4)
  ), nrow = 2, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 2, ncol = 1)
K = both_city[[5]]$nparams
n = both_city[[5]]$nobs
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# CARE: summer vs. winter
R = matrix(
  data = c(
    -1, 1, 0, 0, rep(0, 4)
  ), nrow = 1, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
K = both_city[[5]]$nparams
n = both_city[[5]]$nobs
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Non-CARE: summer vs. winter
R = matrix(
  data = c(
    0, 0, -1, 1, rep(0, 4)
  ), nrow = 1, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Summer: CARE vs. Non-CARE
R = matrix(
  data = c(
    1, 0, -1, 0, rep(0, 4)
  ), nrow = 1, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Winter: CARE vs. Non-CARE
R = matrix(
  data = c(
    0, 1, 0, -1, rep(0, 4)
  ), nrow = 1, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
K = both_city[[5]]$nparams
n = both_city[[5]]$nobs
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Winter: CARE vs. Non-CARE
R = matrix(
  data = c(
    0, 1, 0, -1, rep(0, 4)
  ), nrow = 1, byrow = TRUE)
V = both_city[[5]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
K = both_city[[5]]$nparams
n = both_city[[5]]$nobs
b = both_city[[5]]$coefficients %>% matrix(ncol = 1)

# Summer vs. winter
R = matrix(
  data = c(
    1, -1, 0, 0
  ), nrow = 1, byrow = TRUE)
V = season_city[[3]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
K = season_city[[3]]$nparams
n = season_city[[3]]$nobs
b = season_city[[3]]$coefficients %>% matrix(ncol = 1)

# CARE vs. Non-CARE
R = matrix(
  data = c(
    1, -1, 0, 0
  ), nrow = 1, byrow = TRUE)
V = care_city[[3]]$cov.scaled %>% solve() %>% solve()
q = matrix(0, nrow = 0, ncol = 1)
K = care_city[[3]]$nparams
n = care_city[[3]]$nobs
b = care_city[[3]]$coefficients %>% matrix(ncol = 1)

W = t(R %*% b) %*% solve(R %*% V %*% t(R)) %*% (R %*% b)
# pchisq(q = W, df = nrow(R), lower.tail = FALSE)
J = nrow(R)
pf(q = W / J, df1 = nrow(R), df2 = n - K, lower.tail = FALSE)